A Local Discontinuous Galerkin approximation for the $p$-Navier-Stokes system, Part I: Convergence analysis

In the present paper, we propose a Local Discontinuous Galerkin (LDG) approximation for fully non-homogeneous systems of $p$-Navier-Stokes type. On the basis of the primal formulation, we prove well-posedness, stability (a priori estimates), and weak convergence of the method. To this end, we propose a new DG discretization of the convective term and develop an abstract non-conforming theory of pseudo-monotonicity, which is applied to our problem. We also use our approach to treat the $p$-Stokes problem.

1. Introduction. We consider the numerical approximation of steady systems of p-Navier-Stokes type, i.e., −div S(Dv) + [∇v]v + ∇q = g − div G in Ω , using a Local Discontinuous Galerkin (LDG) scheme. The physical problem motivating this study is the steady motion of a homogeneous, incompressible fluid with sheardependent viscosity. More precisely, for a given vector field g: Ω → R d and a given tensor field G: Ω → R d×d , jointly describing external body forces, as well as a given divergence g: Ω → R, and a given Dirichlet boundary condition v 0 : ∂Ω → R d , we seek for a velocity Mathematical investigation of fluids with shear-dependent viscosities started with the celebrated work of O. Ladyzhenskaya (cf. [44]). In recent years, there has been enormous progress in the understanding of this problem, and we refer the reader to [46,47,7,8] and the references therein for a detailed discussion.
Introducing the unknowns L, S, K : Ω → R d×d , the system (1.1) can be re-written as a "first order" system, i.e., L = ∇v, S = S(L sym ), where we used in (1.3) 2 that [∇v]v = 1 2 div (v ⊗ v) + 1 2 (∇v − div v I d )v and (1.1) 1,2 . Note that this is a modification of the standard Temam modification of the convective term taking into account that the prescribed divergence is non-trivial. Discontinuous Galerkin (DG) methods for elliptic problems have been introduced in the late 90's. They are by now well-understood and rigorously analyzed in the context of linear elliptic problems (cf. [2] for the Poisson problem). In contrast to this, only few papers treat p-Laplace type problems using DG methods (cf. [18,17,30,23,49,52]), or nonconforming methods (cf. [6]). There are several investigations of p-Stokes problems using FE methods (cf. [56,4,3,9,37]) and using DG methods (cf. [19,24,36,13]). Several papers treat the Navier-Stokes problem (p = 2) by means of DG methods. We refer to [21,22,51,25,26,20] to name a few. In addition, the convergence analysis for the p-Navier-Stokes problem using FE methods can be found in [29,43,41,42,34]. On the other hand, to the best of the authors knowledge, there are no investigations using DG methods for the p-Navier-Stokes problem (1.1) and it is desirable to fill this lacuna. This motivated us to develop an abstract framework for the convergence analysis of DG formulations of (1.1) and related problems. This framework extends the quasi non-conforming setting from [10] to the fully non-conforming setting (cf. [39] for the evolutionary case). We exemplify the potential of this abstract framework on an equal order LDG scheme of (1.1) with discontinuous velocity and continuous pressure for low regularity data from (W 1,p 0 (Ω)) * and prove the first convergence for a DG formulation of (1.1). On the other hand, the framework is general enough to be used for other DG schemes, which exploit the advantages of DG methods as e.g. flexibility of the local polynomial degree, flexibility of the meshes, possibility of low order methods, enforcing the divergence constraint exactly with low order elements, or using discontinuous pressure. Moreover, using hybridisation DG methods become competitive in terms of computational costs. None of these aspects is discussed in this paper, but is left for future work.In Part I of the paper, we develop an abstract theory for steady non-conforming pseudo-monotone problems and propose a new DG discretization of the convective term. Moreover, we present easily verifiable sufficient conditions, which ensure the applicability of the abstract theory. Finally, we apply this theory to the p-Navier-Stokes and the p-Stokes problem. For the former, we restrict ourselves to the case p ∈ (2, ∞), since we deal with fully non-homogeneous boundary conditions, while for the latter one, all p ∈ (1, ∞) are covered. In Part II of the paper, we prove convergence rates for the velocity of the homogeneous p-Navier-Stokes problem and the homogeneous p-Stokes problem. In particular, we prove optimal linear convergence rates for the velocity of an LDG scheme with only linear ansatz functions, which is not possible with conforming inf-sup stable FE methods in general shape-regular meshes. Part III of the paper provides convergence rates for the pressure of the same problems. In particular, we prove linear convergence rates for the pressure of the same LDG scheme.
This paper is organized as follows: In Section 2, we introduce the employed notation, define relevant function spaces, basic assumptions on the extra stress tensor S and its consequences, weak formulations in Problem (Q) and Problem (P) of the system (1.1), and discrete operators. In Section 3, we define our numerical fluxes and derive the flux and the primal formulation, i.e, Problem (Q h ) and Problem (P h ), of the system (1.1). In Section 4, we introduce a general concept of non-conforming approximations and nonconforming pseudo-monotonicity. In Section 5, we prove the existence of DG solutions (cf. Proposition 5.6), the stability of the method, i.e., a priori estimates (cf. Proposition 5.7), and the convergence of DG solutions (cf. Theorem 5.8, Theorem 5.29). Theorem 5.8 is the first (weak) convergence result for an LDG method for steady systems of p-Navier-Stokes type. In Section 6, we confirm our theoretical findings via numerical experiments.

Preliminaries.
2.1. Function spaces. We employ c, C > 0 to denote generic constants, that may change from line to line, but are not depending on the crucial quantities. Moreover, we write f ∼ g if and only if there exists constants c, C > 0 such that c f ≤ g ≤ C f .
For k ∈ N and p ∈ [1, ∞], we employ the customary Lebesgue spaces (L p (Ω), · p ) and Sobolev spaces (W k,p (Ω), · k,p ), where Ω ⊆ R d , d ∈ {2, 3}, is a bounded, polyhedral Lipschitz domain. The space W 1,p 0 (Ω) is defined as those functions from W 1,p (Ω) whose trace vanishes on ∂Ω. We equip W 1,p 0 (Ω) with the norm ∇ · p . We do not distinguish between function spaces for scalar, vector-or tensor-valued functions. However, we denote vector-valued functions by boldface letters and tensorvalued functions by capital boldface letters. The standard scalar product between two vectors is denoted by a · b, while the Frobenius scalar product between two tensors is denoted by A : B. The mean value of a locally integrable function f over a measurable set M ⊆ Ω is denoted by |M | M f dx. Moreover, we employ the notation (f, g) := Ω f g dx, whenever the right-hand side is well-defined.
Remark 2.3. (i) Let ϕ be defined in (2.2) and let {ϕ a } a≥0 be the corresponding family of the shifted N-functions. Then, the operators S a : R d×d → R d×d sym , a ≥ 0, defined, for every a ≥ 0 and A ∈ R d×d , via have (p, δ + a)-structure (cf. Remark 2.1). In this case, the characteristics of S a depend only on p ∈ (1, ∞) and are independent of δ ≥ 0 and a ≥ 0.
(ii) Note that the (p, δ)-structure contains as particular cases, e.g., the Ladyzhenskaya model, the Smagorinski model, power law models, the Carreau-Yasuda model and the Powell-Eyring model.

2.3.
The p-Navier-Stokes system. Let us briefly recall some well-known facts about the p-Navier-Stokes system (1.1). For p ∈ (1, ∞), we define the function spaces With this notation, the weak formulation of problem (1.1) is the following.
Problem (Q). For given (g, Alternatively, we can reformulate Problem (Q) "hiding" the pressure. Problem (P). For given (g, where for every f ∈ L p 0 (Ω) The notions "Problem (Q)" and "Problem (P)" are traditional, cf. [16]. Note that In what follows, we always denote by T h , h > 0, a family of uniformly shape regular and conforming triangulations of Ω ⊆ R d , d ∈ {2, 3}, cf. [15], each consisting of d-dimensional simplices K. The parameter h > 0, refers to the maximal mesh-size of T h , i.e., if we define h K := diam(K) for all K ∈ T h , then h := max K∈T h h K . For simplicity, we always assume that h ≤ 1. For a simplex K ∈ T h , we denote by ρ K > 0, the supremum of diameters of inscribed balls. We assume that there is a constant ω 0 > 0, independent of h > 0, such that h K ρ −1 K ≤ ω 0 for every K ∈ T h . The smallest such constant is called the chunkiness of (T h ) h>0 . By Γ i h , we denote the interior faces and put Γ h := Γ i h ∪∂Ω. We assume that each simplex K ∈ T h has at most one face from ∂Ω. We introduce the following scalar products on Γ h if all the integrals are well-defined. Similarly, we define the products ·, · ∂Ω and ·, · Γ i h .

Broken function spaces and projectors.
For every m ∈ N 0 and K ∈ T h , we denote by P m (K), the space of polynomials of degree at most m on K. Then, for given p ∈ (1, ∞) and k ∈ N 0 , we define the spaces (2.10) In addition, for given k ∈ N 0 , we define Q k h,c := Q k h ∩ C 0 (Ω) and V k h,c := V k h ∩ C 0 (Ω). Note that W 1,p (Ω) ⊆ W 1,p (T h ) and V k h ⊆ W 1,p (T h ). We denote by Π k h : Analogously, we define the (local) L 2 -projection into X k h , i.e., Π k h : . For each face γ ∈ Γ h of a given simplex K ∈ T h , we define this interior trace by tr K γ (w h ) ∈ L p (γ). Then, for every w h ∈ W 1,p (T h ) and interior faces γ ∈ Γ i h shared by adjacent elements K − γ , K + γ ∈ T h , we denote by the average and normal jump, resp., of w h on γ. Moreover, for every w h ∈ W 1,p (T h ) and boundary faces γ ∈ ∂Ω, we define boundary averages and boundary jumps, resp., via where n : ∂Ω → S d−1 denotes the unit normal vector field to Ω pointing outward. Analogously, we define {X h } γ and X h n γ for all X h ∈ X k h and γ ∈ Γ h . Furthermore, if there is no danger of confusion, we will omit the index γ ∈ Γ h , in particular, if we interpret jumps and averages as global functions defined on the whole of Γ h .

DG gradient and jump operators.
For every k ∈ N 0 and face γ ∈ Γ h , we define the (local) jump operator R k h,γ : In addition, for every k ∈ N 0 , the DG gradient operator G k h : In particular, for every w h ∈ W 1,p (T h ) and X h ∈ X k h , we have that , the local symmetric gradient. In addition, for every k ∈ N 0 and X k,sym In particular, for every The following discrete Korn inequality on V k h shows that · D,p,h ∼ · ∇,p,h on V k h and, thus, forms a cornerstone of the numerical analysis of the p-Navier-Stokes system (1.1).

Proposition 2.4 (Discrete Korn inequality). For every
Proof. In principle, we proceed as in [14,Lemma 1], where the case p = 2 is treated, and solely replace [14,Inequality (12)] by [13,Proposition 14] in doing so. However, for the convenience of the reader, we give a thorough proof here.
Let N k h be the set of degrees of freedom associated with V k h,c and T h (z) the set of simplices sharing z if z ∈ N k h . For every k ∈ N, the node-averaging quasi-interpolant [30,Lemma A.9]) and an element-wise application of the trace theorem.

From (2.26) and (2.27), using
For the symmetric DG norm, there holds a similar relation like (2.21).
Proof. Appealing to [30, (A.25)], there exists a constant c R > 0 such that for In particular, for every w h ∈ W 1,p (T h ) and z h ∈ Q k h , we have that , and q h := q(q h ), resp., and obtain, also using that For given boundary data v 0 ∈ W 1− 1 p ,p (∂Ω) and given divergence g ∈ L p (Ω) such that (1.2) is satisfied, we denote by v * ∈ W 1,p (Ω) an extension of v 0 with div v * = g in L p (Ω), which exists according to [50, (3.13)]. Using this extension, and restricting ourselves to the case that q h ∈Q k h,c :=Q k h ∩ C 0 (Ω) the numerical fluxes are, for every stabilization |} is defined as in (2.6). Thus we arrive at an inf-sup stable system without using a pressure stabilization. If one wants to work with a discontinuous pressure one has to modify the fluxes as follows: Proceeding as in [30] and, in addition, using that Π k h is self-adjoint, we arrive at the flux formulation of (1.1): For given v * ∈ W 1,p (Ω), g ∈ L p (Ω), G ∈ L p (Ω), and Next, we eliminate in the system (3.7) the variables L h ∈ X k h , S h ∈ X k h and K h ∈ X k h to derive a system only expressed in terms of the two variables v h ∈ V k h and q h ∈Q k h,c .
Note that it follows from (3.7) 1,2,3 that If we insert this into (3.7) 4 , we get the discrete counterpart of Problem (Q): Next, we eliminate in the system (3.8), (3.9), the variable q h ∈Q k h,c to derive a system only expressed in terms of the single variable which is non-empty due to the surjectivity of the Bogovskiȋ operator from W 1,p 0 (Ω) into L p 0 (Ω) and (2.31).
, we obtain the discrete counterpart of Problem (P): (3.10) In particular, note that . Problem (Q h ) and Problem (P h ) are called primal formulations of the system (1.1). In order to formulate Problem (P h ) as an equivalent operator equation in V k h (0), we make the ansatz Using (3.11), the system (3.10) is equivalent to seeking u h ∈ V k h (0) such that for all where, introducing the abbreviations the non-linear operator The derivation of well-posedness (i.e., solvability), stability (i.e., a priori estimates), and weak convergence of the operator equation (3.12) to a weak formulation of the system (1.1) naturally motivates an abstract weak convergence theory based on a non-conforming generalization of the standard notion of pseudo-monotonicity.
4. Non-conforming pseudo-monotonicity. In this section, we introduce a general concept of non-conforming pseudo-monotonicity. It can be understood as a further generalization of the recently introduced quasi-non-conforming pseudo-monotonicity, cf. [10].
Definition 4.1 (Non-conforming approximation). Let V be a reflexive Banach space and (X n ) n∈N a sequence of Banach spaces such that V ⊆ X n with · V = · Xn on V for all n ∈ N and let Y be a Banach space such that X n ⊆ Y for all n ∈ N. Then, a sequence of reflexive Banach spaces V n ⊆ X n , n ∈ N, is called a non-conforming approximation of V with respect to (X n ) n∈N and Y , if the following conditions are satisfied: from sup n∈N v n Xm n < ∞, it follows the existence of a subsequence (v n ) ∈N of (v n ) n∈N and an element v ∈ V such that v n v in Y ( → ∞).
Condition (NC.1) represents the strong approximability of the spaces (V n ) n∈N of V , which can be found in the literature (cf. [32]). Condition (NC.2) is the non-conforming analogue of the weak compactness of bounded sets in reflexive Banach spaces (cf. [26] for a particular realization). The abstract condition seems to be new, since previous convergence results for DG schemes are mostly based on convergence rates, for which this weak condition is not needed.
Next, we introduce an extension of the notion of pseudo-monotonicity to the framework of non-conforming approximations.
Definition 4.2 (Non-conforming pseudo-monotonicity). Let (V n ) n∈N be a nonconforming approximation of V with respect to (X n ) n∈N and Y . Then, a sequence of operators A n : X n → X * n , n ∈ N, is called non-conforming pseudo-monotone with respect to (V n ) n∈N and A :

2)
it follows that Av, v−z V ≤ lim inf n→∞ A mn v n , v n − z Xm n for all z ∈ V. In particular, note that (4.1) already implies that v ∈ V due to (NC.2).
A useful property of pseudo-monotone operators is its stability under summation. This property is shared by non-conforming pseudo-monotone operators. Lemma 4.3. Let (V n ) n∈N be a non-conforming approximation of V with respect to (X n ) n∈N and Y . If A n , B n : X n → X * n , n ∈ N, are non-conforming pseudo-monotone with respect to (V n ) n∈N and A, B : V → V * , resp., then A n + B n : X n → X * n , n ∈ N, is non-conforming pseudo-monotone with respect to (V n ) n∈N and A + B: V → V * .
Proof. Let v n ∈ V mn , n ∈ N, where (m n ) n∈N ⊆ N with m n → ∞ (n → ∞), be a sequence satisfying (4.1) and lim sup n→∞ (A mn + B mn )v n , v n − v Xm n ≤ 0. For n ∈ N, set a n := A mn v n , v n − v Xm n and b n := B mn v n , v n −v Xm n . Then, lim sup n→∞ a n ≤ 0 and lim sup n→∞ b n ≤ 0. In fact, suppose the contrary, e.g., that lim sup n→∞ a n = a > 0. Then, there exists a subsequence such that a n → a ( → ∞) and, consequently, lim sup →∞ b n ≤ lim sup →∞ a n + b n − lim →∞ a n ≤ −a < 0, i.e., a contradiction, as then the non-conforming pseudo-monotonicity of B n : X n → X * n , n ∈ N, choosing z = v, implies 0 ≤ lim inf n→∞ b n < 0. Therefore, lim sup n→∞ a n ≤ 0 and lim sup n→∞ b n ≤ 0 hold, and the non-conforming pseudo-monotonicity of the sequences A n , B n : Summing these inequalities yields the assertion.
The following theorem gives sufficient conditions on a family of non-conforming pseudo-monotone operators such that the associated non-conforming approximation scheme is well-posed, stable, and weakly convergent to a solution of the problem to be approximated.
Theorem 4.4. Let (V n ) n∈N be a non-conforming approximation of V with respect to (X n ) n∈N and Y . Moreover, let A n : X n → X * n , n ∈ N, be non-conforming pseudomonotone with respect to (V n ) n∈N and A : V → V * with the following properties: (AN.1) For every n ∈ N, A n : X n → X * n is pseudo-monotone. (AN.2) There exists a weakly coercive mapping C : R ≥0 → R, i.e., C (a) → ∞ (a → ∞), and a constant c > 0 such that for all z n ∈ V n and n ∈ N, it holds A n z n , z n Xn ≥ C ( z n Xn ) z n Xn − c .
(AN.3) There exists a non-decreasing mapping B : R ≥0 → R ≥0 such that for all z n ∈ V n and n ∈ N, it holds A n z n X * n ≤ B( z n Xn ) .
Let f ∈ V * be a linear functional such that there exists a sequence f n ∈ X * n , n ∈ N, with the following properties: . Then, the following statements apply: (i) For every n ∈ N, there exists v n ∈ V n such that (id Vn ) * A n v n = f n in V * n , i.e., for every z n ∈ V n , it holds A n v n , z n Xn = f n , z n Vn .
Proof. (i). Due to the conditions (AN.1)-(AN.3), for every fixed n ∈ N, the operator A n : X n → X * n is pseudo-monotone, coercive and bounded, which, implies that also the restriction (id Vn ) * A n : V n → V * n is pseudo-monotone, coercive and bounded. As a result, since the spaces (V n ) n∈N are reflexive, the main theorem on pseudo-monotone operators, cf. [57], yields the existence of v n ∈ V n such that (id Vn ) * A n v n = f n in V * n . (ii). For every n ∈ N, using (AN.2) and (i), we find that n v n Xn . Thus, (BN.1) and the weakly coercivity of C : R ≥0 → R imply sup n∈N v n Xn < ∞. Using (AN.3), we conclude that also sup n∈N A n v n X * n < ∞. (iii). Appealing to (ii) and (NC.2), there exists a subsequence v mn ∈ V mn , n ∈ N, where (m n ) n∈N ⊆ N with m n → ∞ (n → ∞), and v ∈ V such that v mn v in Y (n → ∞). In addition, (ii) implies sup n∈N (id V ) * A n v n V * ≤ sup n∈N A n v n X * n < ∞ 2 , so that, by the reflexivity of V * , there exist a not relabeled subsequence and some h ∈ V * such that Let z ∈ D be arbitrary. Then, (NC.1) yields a sequence z n ∈ V n , n ∈ N, such that z n − z Xn → 0 (n → ∞). In addition, (NC.2) yields that z n z in Y (n → ∞), so that, resorting to (BN.2), we deduce that f mn , z mn Xm n → f, z V (n → ∞) . (

4.4)
On the other hand, taking into account (ii) and (4.3), we deduce that (4.5) Since A mn v mn , z mn Xm n = f mn , z mn Xm n for all n ∈ N (cf. (i)) and z ∈ D was chosen arbitrary, combining (4.4) and (4.5), we find that f, z V = h, z V for all z ∈ D, i.e., due to the density of D in V , it holds f = h in V * . Thus, using (BN.2) and (4.3), we get so that the non-conforming pseudo-monotonicity of A n : X n → X * n , n ∈ N, (BN.2), (ii) and (4.3), for every z ∈ V , imply that Eventually, choosing z = v ±z ∈ V for arbitraryz ∈ V in (4.6), we conclude that Av − f,z = 0 for allz ∈ V , i.e., Av = f in V * .
Proof. Let us first verify (BN.1). To this end, let w hn ∈ W 1,p (T hn ) be arbitrary. Then, we have that f k hn , w hn W 1,p (T hn ) ≤ g p w hn p + G p G k hn w hn p . (5.22) Appealing to [30,Lemma A.9], it holds w hn p ≤ c w hn ∇,p,hn , so that from (5.22), using the norm equivalence (2.21), it follows that f k hn (W 1,p (T hn )) * ≤ c g p + c G p .
Now we have prepared everything to prove the well-posedness, stability, and (weak) convergence of our schemes.
Remark 5.9. To the best of the authors knowledge, this is the first convergence results of a DG method for the p-Navier-Stokes problem (1.1). A similar result for the Navier-Stokes problem (p = 2) can be found in [25,26], where a different discretization of the convective term is treated. We would like to emphasize that we treat the fully non-homogeneous problem, which is the reason that we restrict ourselves to the case p > 2. In the homogeneous case, i.e., v 0 = 0 and g = 0, the same methods would work for p > 3d d+2 . Even for the Navier-Stokes problem (p = 2) most previous results (except [25,26]) prove convergence of the method through appropriate convergence rates under additional regularity assumptions on the solution of the limiting problem. In the second part of this paper, we also show convergence rates for our method in the homogeneous case and for G = 0.
This method of proof of course also works for the p-Stokes problem, i.e., we neglect the convective term in (1.1), Problem (P), Problem (P h ), Problem (Q), and Problem (Q h ). For every p ∈ (1, ∞), the existence of solutions for the corresponding discrete problems follows from Lemma 5.2, Lemma 5.5 and Theorem 4.4 (i) in the same way as in the proof of Proposition 5.6. Also the assertions of Proposition 5.7 hold for the p-Stokes problem analogously for every p ∈ (1, ∞), but now with the constant c(v * ) := v * p 1,p . Concerning the convergence of the discrete solutions, we have the following result: Theorem 5.10 (Convergence). Let (h n ) n∈N ⊆ R >0 be such that h n → 0 (n → ∞) and let p ∈ (1, ∞). For every n ∈ N, let (v hn , q hn ) ∈ V k hn ×Q k hn,c be a solution of (3.8), (3.9) without the terms coming from the convective term. Then, there exists a subsequence (h mn ) n∈N ⊆ R >0 , where (m n ) n∈N ⊆ N with m n → ∞ (n → ∞), and (v, q) ∈ V (g) ×Q such that where p * := dp d−p if p < d and p * := ∞ if p ≥ d, and v ∈ V (g) solves Problem (P) without convective term, cf. (2.9), while (v, q) ∈ V (g) ×Q solves Problem (Q) without convective term, cf. (2.7), (2.8).
Proof. We proceed exactly as in the proof of Theorem 5.8 until (5.37), noting that, due to the absence of the convective term, all arguments work for p ∈ (1, ∞). For p > 2, we can finish the proof as before, showing D k hn v hn → Dv in L p (Ω) (n → ∞). For p ≤ 2, we proceed slightly differently. Appealing to Assumption 2.2 and Remark 2.1, it holds This implies, as in [54,Lemma 18], that L sym hn → Dv a.e. in Ω, which together with (2.5), (5.29) and the classical result in [35,Lemma 1.19] yields S(L sym hn ) S(Dv) in L p (Ω) (n → ∞), so that we conclude as in the proof of Theorem 5.8.

Numerical experiments.
In this section, we apply the LDG scheme (3.7) (or (3.8)-(3.9) and (3.10)) to solve numerically the system (1.1) with S : R d×d → R d×d , for every A ∈ R d×d defined by S(A) := (δ + |A sym |) p−2 A sym , where δ := 1e−4 and p > 2. We approximate the discrete solution v h ∈ V k h of the non-linear problem (3.7) deploying the Newton solver from PETSc (version 3.17.3), cf. [45], with an absolute tolerance of τ abs = 1e−8 and a relative tolerance of τ rel = 1e−10. The linear system emerging in each Newton step is solved using a sparse direct solver from MUMPS (version 5.5.0), cf. [1]. For the numerical flux (3.4), we choose the fixed parameter α = 2.5. This choice is in accordance with the choice in [30, Table 1]. In the implementation, the uniqueness of the pressure is enforced via a zero mean condition.
All experiments were carried out using the finite element software package FEniCS (version 2019.1.0), cf. [45].
We construct a initial triangulation T h0 , where h 0 = 1 √ 2 , by subdividing a rectangular cartesian grid into regular triangles with different orientations. Finer triangulations T hi , i = 1, . . . , 5, where h i+1 = hi 2 for all i = 1, . . . , 5, are obtained by regular subdivision of the previous grid: Each triangle is subdivided into four equal triangles by connecting the midpoints of the edges, i.e., applying the red-refinement rule, cf. [5,Definition 4.8 (i)].
Then, for the resulting series of triangulations T hi , i = 1, . . . , 5, we apply the above Newton scheme to compute the corresponding numerical solutions (v i , L i , S i , K i ) := (v hi , L hi , S hi , K hi ) ∈ V k hi × X k hi × X k hi × X k hi , i = 1, . . . , 5, and the error quantities A sym for all A ∈ R d×d . As estimation of the convergence rates, the experimental order of convergence (EOC) where for any i = 1, . . . , 5, we denote by e i either e L,i , e S,i , e ,i , or e q,i , resp., is recorded. For p ∈ {2.2, 2.5, 3, 3.5}, ρ ∈ {0.01, 0.05, 0.1} and a series of triangulations T hi , i = 1, . . . , 5, obtained by regular, global refinement as described above, the EOC is computed and presented in Table 1, Table 2, Table 3 and Table 4, resp. In each case, we observe the expected a convergence ratio of about EOC i (e i ) ≈ ρ p 2 , i = 1, . . . , 5, indicating the convergence of LDG scheme (3.7) (or (3.8), (3.9) and (3.10)). We have also carried out numerical experiments for ρ = 0. Even in this case, the error for all quantities decreases and we see a positive but decreasing very small convergence rate. This indicates that the algorithm converges even in this case, confirming the assertions of Theorem 5.8.  Table 4: Experimental order of convergence: EOCi(eq,i), i = 1, . . . , 5.